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ABSTRACT 


We  are  interested  in  developing  digital  image  restoration 
algorithms  using  a  class  of  spatial  interaction  models  known 
as  conditional  Markov  models.  Our  approach  is  to  represent 
the  images  by  Markov  models  on  toroidal  lattices  and  develop 
minimum  mean  square  error  (MMSE)  restoration  algorithms.  The 
algorithms  are  non-recursive  in  structure  and  due  to  the  under¬ 
lying  representation  on  toroidal  lattices  can  be  implemented 
using  FFT  algorithms. 

We  give  two  types  of  algorithms.  First  we  assume  that  a 
prototype  of  the  original  image  is  available  and  develop 
algorithms  for  restoration  of  degraded  images.  The  degradation 
is  assumed  to  be  due  to  a  known  space  invariant,  non-separable , 
periodic  point  spread  function  and  additive  white  noise. 
Secondly,  we  consider  the  more  general  situation  when  a  proto¬ 
type  image  is  not  available  and  give  algorithms  for  MMSE 
filtering  of  noisy  images.  Experimental  results  are  given  for 
the  above  cases. 
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J .  Introduction 


The  restoration  of  degraded  images  has  many  fields  of 
application,  including  space  imagery  and  biomedical  images 
[1,2].  The  literature  on  image  restoration  is  too  enormous 
to  be  listed  in  detail.  There  are  many  different  methodologies 
such  as  MMSE  restoration,  maximum  a  posteriori  probability 
restoration,  and  maximum  entropy  restoration  [3].  We  are 
concerned  with  MMSE  methods  in  this  paper. 

Suppose  y  and  x  represent  the  lexicographically  ordered 
arrays  of  the  original  and  degraded  image,  related  by 
;;  =  Hy  +  n  (1.1) 

In  (1.1),  H  is  the  block-circulant  matrix  corresponding  to 
the  blur  caused  by  a  non-separable,  space  invariant  point 
spread  function  (PSF)  and  n  is  signal  independent  additive 
white  noise  of  variance  y.  Then  it  is  well  known  [l,p-136] 
that  £,  the  MMSE  estimate  of  y,  can  be  written  as 
?=  Qy  HT(H  Qy  HT  +  YI)_1x  (1.2) 

where  Qy  is  the  covariance  matrix  of  y. 

There  are  two  problems  to  be  considered  in  the  computation 
of  J,  namely,  the  determination  of  the  covariance  matrix  Qy  and 
the  inversion  of  the  matrix  in  (1.2).  Suppose  we  assume  that 
the  prototype  of  the  original  is  available.  Then  Qy  can  be 
evaluated  by  making  appropriate  assumptions  about  the  correlation 
structure  of  the  prototype  image.  For  instance,  one  may  consider 
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a  separable  correlation  function  [ 4 | ,  whose  parameters  are 
chosen  arbitrarily.  In  [1,  Ch.  7],  J  of  (1.2)  is  equivalently 
written  in  terms  of  the  spectral  density  function  (SDF)  of 
the  prototype,  which  is  often  estimated  from  the  periodogram. 

The  second  problem  is  concerned  with  the  matrix  inversion 

2  2 

in  (1.2).  For  an  image  of  size  MxM,  is  of  dimension  M  xM 
and  it  is  not  uncommon  to  have  M  =  128  in  typical  restoration 
applications.  Thus  some  assumptions  have  to  be  made  to  reduce 
the  computational  load.  In  [1,  Ch.  7]  the  block  Toeplitz 
covariance  matrix  is  approximated  by  a  block-circulant 
matrix.  Since  block-circulant  matrices  possess  an  eigen¬ 
function  expansion  in  terms  of  Fourier  vectors,  FFT  computa¬ 
tions  are  used  for  the  implementation  of  (1.2).  In  [4], 
appropriate  assumptions  are  made  regarding  the  model  represen¬ 
tation  for  y,  so  that  has  a  symmetric,  tri-diagonal  Toeplitz 
form.  Since,  symmetric,  tri-diagonal  matrices  are  diagonalized 
by  sine  transforms,  fast  algorithms  for  the  implementation  of 
(1.2)  have  been  developed.  Another  approach  [5,6]  is  to  assume 
an  underlying  unilateral  model  for  y  displayed  as  a  state  space 
model.  Then  y  could  be  computed  by  recursive  algorithms  of 
the  Kalman  type. 

In  this  paper  we  present  a  class  of  MMSE  restoration 
algorithms  which  do  not  make  restrictive  assumptions  such  as 
causality,  isotropy,  or  separable  correlation  function.  We 
assume  that  the  images  are  represented  by  a  class  of  models 


known  as  spatial  interaction  models.  One  of  the  important 
characteristics  of  an  image  is  the  statistical  dependence  of 
the  gray  levels  at  a  lattice  point  on  those  at  its  neighbors. 
Spatial  interaction  models  (familiarly  known  as  random  field 
models)  characterize  this  statistical  dependency  by  representing 
y(i,j),  the  gray  level  at  location  (i,j),  as  a  linear  com¬ 
bination  of  the  gray  levels  at  the  neighboring  pixels  and 
additive  noise.  We  characterize  this  dependency  using  a 
neighbor  set  N  which  is  a  finite  collection  of  pairs  of 
integers  (i^,j.),  not  including  (0,0),  such  that  y ( i , j )  is 
dependent  on  (y(i  +  i^,  j  +  j x)  ,  (ijyj.^  GN}  .  Among  the 
spatial  interaction  models  we  are  primarily  interested  in  the 
class  of  conditional  Markov  (CM)  models  [7-12].  We  represent 
either  the  given  degraded  image,  x  in  (1.1),  or  the  original 
y,  by  appropriate  CM  models  and  formulate  the  MMSE  restoration 
problem.  The  representation  on  toroidal  lattices  leads  to  a 
covariance  matric  Qx  that  is  block-circuiant ,  leading  to  fast 
implementation  of  filters  using  FFT  algorithms. 

Our  contributions  can  be  summarized  as  follows:  First, 


the  CM  models  used  in  this  paper  are  more  general  than  the 
models  used  in  [4].  We  do  not  make  either  separability  or 
isotropy  assumptions.  The  recursive  algorithms  given  in  [5-6] 
are  exact  only  when  the  SDF  of  the  underlying  CM  model  factorizes 


In  general  the  SDFs  of  CM  models  do  not  factorize.  For 
these  models,  the  MMSE  criterion  naturally  leads  to  non¬ 
recursive  algorithms  and  the  recursive  algorithms  are  only 
approximate  since  these  algorithms  use  unilateral  models 
whose  SDF  is  only  an  approximation  to  the  SDF  of  the  true 
CM  model. 

Secondly,  we  give  results  of  filtering  noisy  images 
using  CM  models  without  requiring  the  availability  of  an 
original  image  [9,13]  or  knowledge  of  the  additive  noise 
variance  y,  an  omnipresent  assumption  in  the  literature  on 
MMSE  image  restoration. 

The  organization  of  the  paper  is  as  follows:  in  Section  2, 
we  give  a  brief  description  of  the  representation  of  images 
using  CM  models.  In  Section  3,  we  develop  two  types  of 
restoration  algorithms  and  give  experimental  results, 
discussion  is  presented  in  Section  4. 


Finally, 


( 


2 .  Representation  of  Images  Using  Conditional  Markov  Models 
Suppose  we  are  given  an  image  fy(s),  s  =  (i,j)£ft), 
ft  =  (s  =  (i,j),  lsi,j?iM}  where  y(s)  is  the  gray  level  of  the 
cell  s.  We  partition  the  finite  lattice  ft  into  mutually 
exclusive  and  totally  inclusive  subsets  ft  ,  the  interior  set, 
and  ftfi,  the  boundary  set,  defined  as 

ftj  =  (s  =  (i,j):  seft  and  (s  +  s'‘)£ft  Vs'€N} 


ftfi  =  (s  =  ( i ,  j ) :  s€Q  and  (s  +  s’)£ft  for  at  least  one  s'^N}. 
The  toroidal  lattice  representation  relevant  to  restoration 
of  finite  images  [8-10,12,13],  is  arrived  at  by  assuming  the 
following  representation  for{y(s)}: 

y(s)  =■  f  o  y (s  +  s'  )  +  /v  e (s)  ,  s eftT  (2.1) 

8  {N  S  1 


y(s)  =  E  0  ,  Y,  (s  +  s' )  +  /ve(s),  s€ft 
s'€N  s  1  I 


(2.2) 


where 


Y1(s)  with  s  =  (k,X,)  =  y  [s  +(i,j)]  if  s  +  (i,j)€ft 


and  (i,j)€N 


(2.3) 


=  yl(k  +  i-1)  mod  M  +  1,  (v  +  i  -1)  mod  M  +  1|  if  s  +  ( i ,  j ) 


and  ( i , j ) ON 


(2.4) 


In  (2.1),  N  is  a  symmetric  neighbor  set  in  that  if  s' ON,  then 
-s' ON.  Further  0  =  0  .  The  noise  sequence  e(s)  is  indepen- 

dent  of  y(r),  s  /  r .  In  particular,  E  (e  (s)  |  y  (r) )  «0 ,  r^s  (2.5) 
which  implies 

E(y (s) |y (r) ,  r€ft,  r  ?  s) 


=  z  0S,  y,(s  +  r) 
s'  ON  1 


(2.6) 


Equation  (2.6)  implies  that  t  y ( * )  \  is  Markov  |7]  with  respect 
to  the  symmetric  neighbor  set  N.  The  noise  sequence  {e(s)l  is 
correlated  with  the  correlation  structure  in  (2.7),  which 
follows  from  (2.5)  or  (2.6): 

E(e(s)  e(s') )  =  1  s  =  s' 

=  -0S_S.  /  s  -  s'€N  (2.7) 

=  0  otherwise 

To  ensure  stationarity ,  (2.8)  should  be  valid  [9]: 

Ms  =  d-2  0T^s)  >  0  (2.8) 

where 

=  Col.  [cos  ((s-l)Ts*  ]  ,  s'  6Ns] 

Here  Ng  is  a  proper  subset  of  N  such  that 

s€n  »  -sj£N 
s  s 

N  UN  =  N 

S  I 

and  N  =  f  s'  :  -s'  €N  } 

—  s 

s 

Using  lexicographically  ordered  notation  for  the  arrays  {y(s)} 
and  (e(s)},  (2.1)  -  (2.2)  can  be  written  in  matrix/vector 
format  as 

A(6)y=/ve  (2.9) 

where  A ( 0 )  is  a  symmetric  block-circulant  matrix.  As  A(0)  has 
an  inverse  by  (2.8), 

y  =  /vA(0) -1e  (2.10) 

and  the  covariance  matrix 

Qy  =  E(yyT)  =  va“1(0)  (2.ii) 

is  also  block-circulant  with  eigenvalues  v/y  ,  s€ft.  By  different 

s 


choices  of  N  we  account  for  different  filter  structures. 
For  instance  N  =  {  (0 , 1)  ,  ( 1 , 0)  ,  (-1 , 0)  ,  (0  ,-D  )  corresponds  to 
the  first  order  Markov  model  [7]. 


In  developing  restoration  algorithms  we  consider  the 


following  two  cases. 

Case  (i) :  A  CM  model  with  parameters  (0,v)  is  assumed  for 
the  original  image  y.  It  is  assumed  that  a  prototype  of 
the  original  image  is  available  such  that  (&,v)  can  be  esti¬ 
mated.  A  degraded  image  x  related  to  y  as  in  (1.1)  is  given 
and  it  is  required  to  develop  algorithms  for  computing  y. 

Case  (ii) ;  We  are  given  a  noisy  image  only,  with  no  prototype 
of  the  original  being  available.  It  is  required  to  develop 
filtering  algorithms. 

3.1. a  Restoration  Using  CM  Models  for  y 
(Case  i) . 

Assume  that  the  original  image  y  is  represented  by  (2.1) 
and  (2.2)  and  that  x  is  related  to  y  as  in  (1.1).  Then  the  MMSE 
estimate  y  is  given  in  (1.2). 


Let  f  =  Col. 

[  t .  ,  A  .  t .  ,  .  . 

.,A  .^t.]. 

with  s  = 

~s 

-D  l  3 

i  ~3 

t  =  Col. 

[1,A.,A.2,.. 

\  M-l.  M 
.  ,  A  .  ]  ,  M 

vector 

~s 

3  3 

3 

A^  =  exp  [/-1A  (i-1)  ]  ,  Ag  =  27i/M 

2 

be  the  M  Fourier  vectors  and  let  h  ,  s€&  be  the  eigenvalues 
of  H.  Then  due  to  the  underlying  toroidal  lattice  representa¬ 
tion,  y  can  be  computed  as 

i  -  iA,£?s('’Ks/MIKsl|2+Y  Ms)ffx  (3.1) 

where  *  denotes  the  complex  conjugate  operator.  The  computation 


2 

complexity  is  0(M  logM) .  The  computing  scheme  in  (3.1)  has 
a  general  structure  without  involving  any  assumptions  such  as 
isotropy  or  separable  correlation  structure.  Equation  (3.1) 
is  characterized  by  v,  0  and  y.  It  is  assumed  that  y  is  known. 
The  parameters  0,v  are  assumed  arbitrarily  [4]  or  estimated 
from  the  prototype  of  the  original.  There  are  several  methods 
of  estimating  (0,v)  for  CM  models.  The  coding  method  in  [11] 
is  very  simple  in  principle  and  yields  consistent  but  very 
inefficient  estimates.  The  ML  method  in  [12]  yields  asympto¬ 
tically  consistent  and  efficient  estimates  but  is  computationally 
unattractive  except  for  isotropic  models.  A  simple  consistent 
estimate  whose  efficiency  lies  between  the  coding  and  ML  methods, 
given  in  [14],  is  used  to  estimate  0  and  v. 


A  numerical  figure  of  merit  for  restoration  algorithms  is 
the  expected  minimum  error  u) .  From  [l,pp.  140], 


a)  =  -„  Tr  [Q  -  L  H  Q  1 


(3.2) 


M 


where 


(3.3) 


L  =  QyHT(H  Qy  ht  +  yi)"1 
Using  the  fact  that  Qy  and  H  are  block-circulants ,  with 


eigenvalues  v/u  and  h  ,  s£ft, 
s  s 


0)  = 

M^ 


E  VY 

ft  v  |E  P" 


(3.4) 


+  YUg) 


3.1.b  Experimental  Results;  To  illustrate  the  use  of 
(3.1),  the  following  experiment  was  done.  A  64x64  window  of 
the  original  USC  girl  image  was  used.  The  parameters  corres¬ 


ponding  to  different  neighbor  sets  N  were  computed.  The 

s 


(3.5) 


image  was  blurred  by  using  the  PSF 

h  (k,  £)  =  5-ii  exp  f  -  0 . 4  (k2  +  £2)},  0  s  k,  £  s  2, 
and  Gaussian  noise  of  signal  to  noise  ratio  (defined  as  the 
ratio  of  signal  variance  to  noise  variance)  of  5  was  added  to 
the  blurred  image.  The  results  of  implementing  (3.1)  are  in 
Fig.  1  for  SNR  =  5  and  in  Fig.  2  for  SNR  =  1.  The  details  of 
the  models  used  are  in  the  second  column  of  Table  I.  The 
restored  images  corresponding  to  neighbor  sets  are  good 

both  for  the  SNR  =  5  and  1  cases.  But  the  restored  images 
corresponding  to  neighbor  sets  Ng^,  N  2  are  not  as  good  for 
SNR  =  1  as  they  are  for  SNR  =  5  blurred  images.  Filters  of 
relatively  large  order  were  used  for  this  image  since  lower 
order  CM  models  were  non-stationary . 

To  get  some  idea  as  to  the  amount  of  theoretical  improve¬ 
ment  that  is  possible  when  using  different  neighbor  sets, 


numerical  values  of  the  expected  error  <o  were  computed  using 


(3.4).  The  quantity  g  defined  as 


g  =  lOlog 


MSE  between  original  and  degraded 


(3.6) 


is  tabulated  in  Table  I  for  SNR  =  5  and  1  for  the  blur  in  (3.5)  . 


The  values  of  the  numerator  of  (3.6)  are  511.41  and  1910.1 


respectively. 

Note  that  thus  far  we  have  assumed  that  the  additive  noise 


sequence  n  is  white.  The  case  when  g  is  colored  [4]  can  be 

•  T 

considered  by  replacing  yl  in  (1.2)  by  C  where  C  =  Ei  g  )  and 


is  known. 


a 


3. 2. a  Restoration  using  CM  models  for  x:  In  practice, 

the  CM  models  of  the  original  image  are  not  knov  \  and  only 

the  degraded  image  is  available.  It  is  required  to  obtain  y 

using  only  the  given  degraded  image.  Suppose  that  degradation 

is  due  to  additive  noise  only,  i.e.,  x  is  related  to  y  by 

x  =  y  +  n  (3.7) 

Assume  that  x  is  represented  by  the  CM  model 

x(s)  =  T  0'  x(s  +  s')  +A?e(s),  s€0  (3.8) 

s  *  I 

and 

x,8)  *  +  •'*>  '*‘aB  (3'9) 

where  x^ ( • )  obeys  an  equation  similar  to  y^(*)  in  (2.3)  and 
(2.4).  The  covariance  matrix  Q  of  x  under  the  transformation 

^  X  -- 

in  (3.7)  is 


0 

=  0 

+  Y I 

(3.10) 

-  X 

~y 

0 

■y 

=  Qx 

-  Y I 

(3.11) 

By  substitution  of  (3.11)  in  (1.2)  and  letting  H  =  I ,  we 
obtain  the  MMSE  estimate 

y  =  (0  ~  Y I )  (Q v)  _1x  (3.12) 

-  A  ~  A  -v 

A 

Since  Q  is  block-circulant  with  eigenvalues  v'/V  ,s€fi,  y  can 

A  S  ~ 


be  computed  as 


•i  (V  -  yy’)  *T 

v  =  £  T.  f  _ f  x 

%  M2  (2  is  v'  Es  X 


(3.13) 


Equation  (3.13)  yields  an  optimal  y  if  the  parameters  6',  v'  and  y 
are  known  exactly.  In  practice  they  are  not  known  and  are 
estimated  from  the  given  noisy  image.  As  discussed  before 
0',  v*'  can  be  estimated  using  the  estimation  scheme  in  [14]  and 


the  noisy  image.  A  scheme  which  has  been  found  to  be  of 
some  practical  use  [1,  pp.  144]  in  estimating  y  is  to  use 
the  steady  state  component  of  the  SDF  of  the  noisy  image  as 
an  estimate.  An  estimate  y  of  y  is  obtained  for  64 x  64  images 
using 

Y  =  jg  I  Sx(s>'  Qu  *  ts  =  (i , j J ,30  s  i , j  s,  33}  (3.14) 

% 

where  s  (s)  =  v’/u',  s is  the  two-dimensional  discrete 

S 

SDF  of  x  under  the  representation  (3.8)  -  (3.9). 

3.2.b  Experimental  Results:  The  original  girl's  face 
in  Fig.  1  was  corrupted  by  additive  Gaussian  noise  of  SNR  =  5. 

It  was  assumed  that  only  the  noisy  image  was  available.  CM 
models  having  different  neighbor  sets  N  were  fitted  to  the 

b 

noisy  image.  The  filtered  images  for  the  SNR  =  5  noisy  image 

obtained  using  (3.13)  are  shown  in  Fig.  3.  The  details  of 

the  models  used  are  given  in  Table  II.  Similar  results  for 

SNR  =  1  are  given  in  Fig.  4.  The  filtered  images  of  the  SNR  =  5 

noisy  image  corresponding  to  neighbor  sets  N  . ,  N  c  and  N  are 

s i  Sj  so 

good.  The  improvement  is  noticeable  in  the  filtered  images  of 
the  SNR  =  1  noisy  image  corresponding  to  neighbor  set  N 

S  X 

The  parameter  y  was  estimated  using  the  scheme  discussed 
in  Section  3. 2. a.  The  best  estimate  of  y  was  obtained  for  the 
SNR  =  5  image;  it  was  378.74  (true  value  being  393.01)  using 
neighbor  set  Ng^,  and  for  a  0  db  noisy  image  the  estimate  was 
1917.9  (the  true  value  being  1965.05)  using  neighbor  set  Ngg. 


Actual  values  of  MSE  between  the  original  and  the  filtered 
images  in  Fig.  3  and  4  were  computed.  The  MSE  corresponding 
to  the  best  filtered  image  in  Fig.  3  was  91.08  using  the 
neighbor  set  Ng^  with  the  error  between  the  original  and  noisy 
image  being  378.53,  roughly  a  reduction  in  error  by  a  factor 
of  4.  Similar  test  results  for  the  SNR  =  1  image  are  321.79 
and  1728.5  using  Ngg,  representing  a  reduction  in  noise  by  a 
factor  of  5.  For  the  details  of  estimates  of  Y  and  MSE  for 
other  CM  models  see  Tables  II  and  III. 


4. 


Discussion 


We  have  considered  some  MMSE  algorithms  for  the 
restoration  of  degraded  images  using  CM  models.  Assuming 
that  a  prototype  of  the  original  image  is  available,  an 
algorithm  for  restoring  a  blurred  and  noisy  image  has  been 
given  together  with  an  expression  for  the  expected  minimum 
error.  A  computing  scheme  for  the  filtering  of  noisy 
images,  when  the  prototype  of  the  original  image  is  not 
available,  has  been  given  along  with  experimental  results. 

The  algorithms  are  fast,  non-recursive,  and  are  valid 

for  any  arbitrary  neighbor  sets  Ng.  We  have  used  a  toroidal 

lattice  representation  for  the  images,  which  might  appear 

restrictive.  For  the  values  of  M  considered  here  it  can  be 

shown  [9,15]  that  the  second-order  properties  of  a  CM  model 

on  a  toroidal  lattice  are  very  nearly  equal  to  those  of  the 

corresponding  CM  model  on  an  infinite  lattice.  To  illustrate 

this  the  normalized  autocorrelation  function  (ACF)  of  an 

isotropic  CM  model  with  8=  .24  and  N  =  {  ( 1 , 0)  ,  (0 , 1)  }  was 

s 

evaluated  at  a  few  lower  lags  using  the  toroidal  assumption 
and  compared  with  the  corresponding  quantities  for  the  infinite 
lattice  model  available  in  [16,  Table  1],  Denoting  by  ^ 
the  normalized  ACF  at  lag  (k,?,)  we  have,  for  the  ACF  computed 
with  infinite  lattice  model,  Q  =  .434,  p.^  ^  =  .295,  PQ  2  =  *220 

and  P2  1  =  -180.  The  corresponding  quantities  computed  for  the 


toroidal  lattice  are  p^  ^  =  .4329,  p^  ^  =  .2935, 
and  p0  .  =  .1785  which  are  numerically  close. 

2.  1  L 
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Table  I .  Performance  bounds  in  decibels  as  predicted 

by  various  CM  models  for  the  PSF  in  (3.5)  and 
M  =  64. 


SNR  =  SNR  = 

Number  Neighbor  set  Ng  7  db  0  db 


{  (1,0)  , (0,1) , (1,1) , (1,-1)  , 

N  5.72  7.44 

Sl  (0,2), (2,0)} 


{  (0,1)  ,  (1,0)  ,  (1,-D  ,  (1,1)  ,  (0,2)  , 

N  5.27  6.66 

s2  (2,0)  ,  (-2,1) , (2,1) , (-1,2) , (1,2) } 


N 


s3 


{  d,0)  ,  (0,1)  ,  (1,-1)  ,  (1,1)  ,  (0,2)  ,  (2,0)  , 
(-2,1)  ,  (2,1)  ,  (-1,2)  ,  (1,2)  ,  (2,2)  ,  (-2,2)  } 


Table  II.  The  estimates  of  y  obtained  using  (3.14).  The 


true  values  of  y  are 
(SNR  =  1) . 

393.01  (SNR  =5) 

and  1965 

Neighbor  set 

SNR  =  5 

SNR  =  1 

s4  =  { (0,1) , (1,0) } 

264.63 

1183.9 

lg5  =  { (0,1) , (1,0) , (1,1) } 

318.62 

1493.0 

ls6  =  ( (0,1) , (1,0) ,  (1,1)  ,  (1,-1)  } 

358.79 

1917.9 

N 


sl 


378.74 


2185.1 


•b* 

& 

_ 4 

Fig.  l.  Restoration 
of  image  containing 
blur  (3*5)  and  Gaussian 
noise  (SNR=5) . 

a)  Original 

b)  Blurred 

c-e)  Restored  images 
using  ,  N  ,  N  . 


Fig.  2.  Analogous  to 
Fig.  1  with  SNR=1. 


Fig.  3.  Filtering  of 
noisy  image  when  proto¬ 
type  of  original  is  not 
available. 

a)  Noisy  image  (SNR=5) . 
b-e)  Filtered  images 
using  Ns4,  Ns5,  N>r 
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